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Abstract 

At  the  mesoscopic  scale,  chemical  processes  have  probability  distributions  that 
evolve  according  to  an  infinite  set  of  linear  ordinary  differential  equations  known 
as  the  chemical  master  equation  (CME).  It  is  commonly  believed  that  the  CME 
cannot  be  solved  except  for  the  most  trivial  of  cases,  but  recent  work  has  raised 
questions  regarding  validity  of  this  belief.  Eor  many  cases,  Einite  State  Projection 
(ESP)  techniques  [1]  can  reduce  the  order  of  the  CME  to  a  solvable  system  while 
retaining  any  prespecified  error  tolerance.  Even  when  accuracy  demands  require  a 
projection  that  is  too  large  to  be  solve  efficiently,  the  ESP  retains  the  linearity  of  the 
CME,  and  is  open  to  a  host  of  additional  model  reductions  and  computational  tech¬ 
niques.  In  this  paper,  we  develop  a  new  algorithm  based  upon  the  linearity  property 
of  super-positioning,  and  we  illustrate  the  benefits  of  this  algorithm  on  a  simplified 
model  of  the  heat  shock  mechanism  in  E.  eoli.  The  new  algorithm  retains  the  full 
accuracy  of  the  original  ESP  algorithm,  but  with  significantly  increased  efficiency 
and  a  greater  range  of  applicability. 
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1  Introduction 


When  studying  processes  at  the  mesoscopic  level,  researchers  often  assume 
that  they  evolve  according  to  a  continuous  time,  jump  Markov  process.  In 
this  regime,  the  system  is  best  described  not  in  terms  of  individual  trajectories 
but  in  terms  of  how  the  system’s  probability  distribution  evolves.  For  chemical 
reactions,  this  distribution  has  been  shown  to  evolve  according  to  the  inhnite 
dimensional  linear  ordinary  differential  equation  (ODE)  known  as  the  chemical 
master  equation  (CME)  [2]. 

It  is  commonly  believed  that  master  equations,  such  as  the  CME,  cannot  be 
solved  except  in  the  simplest  of  circumstances.  This  belief  has  driven  the  com¬ 
putational  research  community  to  devise  clever  kinetic  Monte  Carlo  (KMC) 
algorithms  to  simulate  system  dynamics.  Gilespie’s  stochastic  simulation  al¬ 
gorithm  (SSA)  is  one  such  algorithm  developed  specihcally  for  the  CME  [3]. 
In  the  SSA  each  reaction  is  simulated  using  two  pseudo-random  numbers-the 
hrst  tells  when  the  next  reaction  will  occur,  and  the  second  determines  which 
reaction  it  will  be.  The  SSA  is  considered  exact  in  the  sense  that  if  one  were 
to  conduct  the  simulation  an  inhnite  number  of  times  with  an  ideal  random 
number  generator,  the  collected  statistical  data  would  converge  to  the  exact 
solution  to  the  CME.  Unfortunately,  the  convergence  rate  for  any  KMC  rou¬ 
tine  is  very  slow;  cutting  the  error  in  half  requires  four  times  the  number  of 
simulations.  Since  a  single  simulation  may  contain  huge  numbers  of  individ¬ 
ual  reactions,  the  SSA  may  take  far  too  long  to  be  of  any  use.  Researchers 
have  greatly  improved  efficiency  of  the  SSA  through  various  approximation 
schemes.  Some  of  these  approximations  are  made  by  separating  the  fast  dy¬ 
namics  from  the  slow  [4-8].  During  short  periods  of  time,  the  fast  dynamics 
dominate,  and  the  slow  dynamics  may  be  ignored.  Eor  long  periods  of  time, 
one  averages  out  the  fast  dynamics  equilibrate  in  order  to  emphasize  the  slow 
dynamics.  Other  approximations  are  made  by  discretizing  the  time  interval 
into  r  leaps  and  approximating  the  dynamics  over  those  subintervals  [9-14]. 
Both  approximation  types,  system  partitioning  and  r  leaping,  have  been  very 
successful  in  increasing  the  scope  of  problems  to  which  KMC  schemes  such  as 
the  SSA  may  be  reasonably  applied. 

We  recently  developed  the  Finite  State  Projection  algorithm  as  an  additional 
tool  for  the  analysis  of  jump  Markov  processes  [1].  Unlike  KMC  algorithms, 
the  FSP  algorithm  solves  the  CME  to  any  prespecihed  degree  of  accuracy 
without  random  number  generation.  The  FSP  approach  is  based  upon  linear 
systems  theory  and  works  by  projecting  the  intractable  inhnite  dimensional 
master  equation  onto  a  solvable  hnite  dimensional  space.  Previous  implemen¬ 
tations  of  the  FSP  method  have  been  very  successful  for  some  biologically 
inspired  systems  [1,15],  but  for  many  problems  the  projection  remains  too 
large  to  solve  efficiently,  and  further  model  reductions  are  needed.  While  some 
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of  these  reductions  come  from  fields  such  as  modern  controls  theory  [16] ,  oth¬ 
ers  are  inspired  from  similar  reductions  applied  to  the  KMC  algorithms.  For 
example,  as  we  show  in  [17,18],  we  can  reduce  the  FSP  significantly  using  time- 
scale  separation  based  system  partitioning  approaches.  In  the  current  work, 
we  further  advance  the  FSP  by  exploring  how  some  tools  of  time  discretization 
can  also  be  used  for  the  dramatic  improvement  of  the  FSP. 

In  the  next  section  we  provide  some  background  on  the  basic  FSP  method.  In 
section  3  we  illustrate  how  linearity  of  the  FSP  allows  us  to  efficiently  deal  with 
initial  probability  density  vectors  that  contain  many  non-zero  elements.  In 
section  4  we  provide  a  multiple  time  step  version  of  the  FSP  and  demonstrate 
its  use  with  an  example  from  the  held  of  systems  biology.  Finally,  in  section 
5,  we  conclude  with  remarks  on  the  advantages  of  these  approaches  over  the 
original  FSP  and  outline  a  few  directions  for  future  work  on  this  topic. 


2  The  Finite  State  Projection  Method 


Although  the  hnite  state  projection  (FSP)  method  presented  here  is  valid 
for  any  continuous  time,  discrete  space  Markov  process,  we  present  it  in  the 
context  of  the  chemical  master  equation. 

We  consider  a  system  of  N  chemically  reacting  species.  The  non-negative 
populations  of  the  N  molecular  species  jointly  dehne  a  unique  conhguration 
of  the  system,  x  ;=  [^j^  ^2  ...  ]^-  ^7  fixing  a  sequence  xi,  X2, . . .  of  elements 

in  we  can  define  the  ordered  configuration  set  as  X  ;=  [xi,  X2, ...  ]^.  Let 

Pi{t)  denote  the  probability  that  the  system  will  have  the  configuration  Xj 
at  time  t.  Under  the  assumptions  that  the  system  is  continually  well-mixed 
and  kept  at  constant  volume  and  temperature,  one  can  show  that  the  system 
evolves  according  to  a  discrete  state,  continuous  time  jump  Markov  process, 
whose  distribution,  P  =  [pi,p2,  ■  ■  •]^,  evolves  according  to  the  linear  ordinary 
differential  equation  known  as  the  chemical  master  equation  (CME)  [19]: 


P(«)=AP(«),  (1) 

The  infinitesimal  generator  matrix.  A,  is  defined  by  the  reaction  stoichiome¬ 
tries  and  propensities  and  the  choice  of  the  enumeration  of  X.  Like  any  gen¬ 
erator  matrix,  A  has  the  properties  that  all  diagonal  elements  of  are  non¬ 
positive;  all  off-diagonal  elements  are  non-negative;  and  all  columns  sum  to 
zero.  When  the  set  X  is  finite  dimensional,  one  can  easily  compute  the  solu¬ 
tion  P(t/)  =  exp(Atj)P(0),  but  when  the  set  X  is  infinite  or  extremely  large, 
the  corresponding  solution  is  unclear  or  vastly  difficult  to  compute.  For  these 
cases,  we  devised  the  Finite  State  Projection  (FSP)  method  [1]. 
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To  best  review  the  FSP  method  and  present  our  current  work,  we  must  first 
introduce  some  convenient  notation.  Let  J  =  {ji,  js,  •  •  •}  denote  an  index 
set  to  which  the  usual  operations  (U,  fl;  etc.)  and  relations  (c,  C,  etc.)  ap¬ 
ply.  Let  J'  denote  the  complement  of  the  set  J.  If  X  is  an  enumerated  set 
{xi,X2,X3, . . then  Xj  denotes  the  subset  {xj^,  x^j,  x^g, . . .}.  Furthermore, 
let  vj  denote  the  subvector  of  v  whose  elements  are  chosen  according  to  J, 
and  let  Ajj  denote  the  submatrix  of  A  such  that  the  rows  have  been  chosen 
according  to  I  and  the  columns  have  been  chosen  according  to  J.  For  example, 
if  /  and  J  are  dehned  as  {3, 1,2}  and  {1,3},  respectively,  then: 


a  h  c 

9  k 

d  e  f 

= 

a  c 

g  h  k 

ij 

df  _ 

For  convenience,  let  Aj  ;=  Ajj.  We  will  also  use  an  embedding  operator, 
Vj{.}  as  follows.  Given  any  vector,  v  and  its  J  indexed  subvector,  vj,  the 
vector  Vj  {vj}  has  the  same  dimension  as  v  and  its  only  non-zero  entries  are 
the  elements  of  vj  distributed  according  to  the  indexing  set  J .  Finally,  we  will 
use  the  vector  e*  to  denote  a  vector  with  a  one  in  the  location  and  zeros 
elsewhere.  With  this  notation  we  can  restate  the  following  two  theorems  from 
[1]: 

Theorem  1  If  A  E  has  no  negative  off-diagonal  elements,  then  for  any 
index  set,  J , 

[expA]j  >  exp{Aj)  >  0.  (2) 


Theorem  2  Consider  any  Markov  proeess  in  whieh  the  distribution  evolves 
aeeording  to  the  linear,  time-invariant  ODE: 

P(«)  =  AP(«), 


If  for  some  finite  index  set  J,  e  >  0,  and  tj  >  0, 

||exp(Ajt/)P, 7(0)11^  >  1  - 


then 

and 


Vj{exp{Ajtf)Pj{0)}  <  P(t/), 
P{tf) -Vj{exp{Ajtf)Pj{0)}\\-^  <  e. 


(3) 

(4) 

(5) 


Using  the  current  notation,  these  theorems  have  been  strengthened  slightly 
from  their  original  form,  but  their  proofs  are  easily  found  with  slight  modifi- 
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cation  to  those  presented  in  [1],  Theorem  1  guarantees  that  as  we  add  points 
to  the  hnite  conhguration  subset,  the  approximate  solution  monotonically  in¬ 
creases,  and  Theorem  2  provides  a  certihcate  of  how  close  the  approximation 
is  to  the  true  solution.  Together  the  two  theorems  suggest  the  FSP  algorithm 
presented  in  Ref  [1]: 

The  Original  Finite  State  Projection  Algorithm 


Inputs  Propensity  functions  and  stoichiometry  for  all  reactions. 
Initial  probability  density  vector,  P(0). 

Final  time  of  interest,  tf. 

Total  amount  of  acceptable  error,  e  >  0. 

Step  0  Choose  an  initial  hnite  set  of  states,  Xj^,  for  the  FSP. 
Initialize  a  counter,  i  =  0. 


Step  1  Use  propensity  functions  and  stoichiometry  to  form  Aj.. 

Compute  Pj.  =  ||exp(Aj.t/)Pj.(0)||^. 

Step  2  If  Pj.  >  1  —  e.  Stop. 

Vj.  {exp(A jij)Pj,(0)}  approximates  P(t/)  to  within  a  total  error  of  e. 


Step  3  Add  more  states  to  hnd  X . 

Increment  i  and  return  to  Step  1. 


For  the  basic  FSP  algorithm,  if  we  wish  to  hnd  a  solution  that  is  accurate 
to  within  e  at  a  time  tf,  we  must  hnd  a  hnite  set  of  conhgurations  such  that 
the  probability  of  ever  leaving  that  set  during  the  time  interval  [0,f/]  is  less 
than  e.  For  many  problems,  including  the  examples  shown  in  [1,15],  this  set  of 
conhgurations  may  be  small  enough  that  we  can  easily  compute  a  single  ma¬ 
trix  exponential  to  approximate  the  solution  to  the  CME.  However,  in  other 
situations  the  conhguration  space  required  for  a  one  matrix  solution  may  be 
exorbitantly  large.  In  this  work  we  utilize  the  linearity  and  time  invariance  of 
the  chemical  master  equation  to  address  two  such  situations.  First,  in  section 
3  we  extend  the  FSP  to  handle  problems  in  which  the  initial  probability  dis¬ 
tribution  is  supported  over  a  large  portion  of  the  conhguration  space.  Second, 
in  section  4  we  address  the  problem  that  occurs  when  non-equilibrium  distri¬ 
butions  tend  to  drift  over  time  and  cover  large  portions  of  the  conhguration 
space. 
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3  The  FSP  for  Non-Sparse  Initial  Distributions 


Although  the  original  FSP  method  is  valid  for  any  initial  probability  distri¬ 
bution,  all  examples  in  previous  work  [1,15-18]  begin  with  a  specific  known 
initial  configuration;  if  the  system  begins  in  state  k,  the  initial  probability 
distribution  for  the  CME  was  written,  Pi(0)  =  5*^,  where  6ik  is  the  Kronecker 
delta.  Suppose  now  that  the  initial  distribution  is  given  not  by  the  Kronecker 
delta  but  by  a  vector  with  many  non-zero  elements.  For  example,  suppose  that 
the  initial  distribution  is  specified  by  the  solution  at  the  end  of  a  previous  time 
step.  From  Theorem  2,  in  order  for  the  original  FSP  algorithm  to  converge,  we 
must  be  able  to  find  a  set  of  states,  Xj,  that  satisfies  the  stopping  criterion: 

||exp(Ajt/)Pj(0)||^  >  (1  -  £)■  (6) 

Since  the  sum  of  the  FSP  solution  at  cannot  exceed  the  sum  of  the  truncated 
initial  pdv,  Pj(0),  we  must  always  include  at  least  as  many  states  in  the  FSP 
solution  as  is  required  such  that  ||Pj(0)||^  >  1  —  £•  For  a  sparse  pdv,  such  as 
that  generated  by  <5*^,  this  restriction  on  the  size  of  the  FSP  solution  is  trivial:  J 
need  only  include  k.  However,  when  the  initial  pdv  has  broad  support,  the  size 
of  the  FSP  solution  may  be  much  larger  and  therefore  require  the  inefficient 
calculation  of  very  high-dimensional  matrix  exponentials.  Fortunately,  one  can 
use  the  property  of  superpositioning  guaranteed  by  the  linearity  of  the  FSP 
to  mitigate  this  concern  and  recover  some  computational  efficiency  as  we  can 
show  in  the  following  proposition. 

Proposition  3  Superposition  of  FSP  Solutions 

Consider  any  Markov  process  in  which  the  distribution  evolves  according  to 
the  linear  ODE: 

Pit)  =  APit). 

Let  7  <  1,  (5  <  1  and  tf  >0.  If  there  is  an  index  set  I  such  that: 


l|i’;(0)lli>7,  (7) 

and  if  for  every  i  E  I,  there  is  a  corresponding  index  set  Ji  containing  i  such 
that 

then, 


expiAj^tf)e\\\^  >  S, 
Y^PiVj^  {exp{Ajdf)e\}  <  P{tf), 


iei 


and 


Pitf)  - 


i&I 


<1  —  7(5. 


(8) 

(9) 

(10) 
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Proof.  We  begin  by  proving  (9).  If  we  define  the  index  set  If  =  Uie/  then 
we  have  the  relation, 

Vi^  {exp(A/^t/)P7^(0)}  =  ^  Pi{0)Vi^  {exp(A7^t/)e^J  ,  (11) 

idif 

Since  I  C  If ,  we  are  guaranteed  that 

Vi^  {exp(A/^t/)P/^(0)}  >  {exp(A/^t/)e^J  .  (12) 

i&I 

Furthermore,  since  for  every  i,  Jj  C  If  and  Pi(0)  >  0,  Theorem  1  guarantees 
that. 


Vi^  {exp(A/^t/)P/^(0)}  >  J2pi{0)Vj^  {exp(Ajij)e*jJ  .  (13) 

i&I 

Furthermore,  using  the  result  from  Theorem  1  that  exp(Ajtj)  is  non-negative 
for  any  index  set  J,  and  applying  conditions  (7)  and  (8)  yields 


Vi^  {exp(A/^tj)Pj^(0)}||^ 


> 


^Pi(0)T>j,  [exp{Ajif)e\} 
iei  1 


><5||P/(0)lk 

>  57. 


(14) 


Theorem  2  tells  us  that 

Vi^  {exp(A/^tj)Pj^(0)}  <  P(t/), 

and  then  from  Eqn  (13)  we  show  that 

^Pi(0)T>j,  [exp{AjIf)e\}  <  P{tf), 
ie/o 

which  is  Eqn.  (9). 

Combining  the  fact  that  ||P(t/)||^  =  1  and  inequality  (14)  gives: 


{exp(Aj,t/)e*jJ  >  (||P(t/)||^  -  1)  +57 

iei  1 


Rearanging  this  result  and  applying  15  yields  inequality  (10) 


P(^/)  -  {exp(Ajij)e*jJ 


i&I 


<  1  —  57 


(15) 


(16) 


(17) 


and  completes  the  proof. 
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The  result  of  Proposition  3  now  enables  us  to  modify  the  above  FSP  algo¬ 
rithm  to  better  handle  situations  in  which  the  initial  probability  distribution 
is  non-sparse.  Before  stating  this  new  algorithm,  however,  we  would  first  like 
to  make  a  few  notes  to  explain  our  choice  of  notation.  First,  although  this 
algorithm  can  be  useful  on  its  own,  we  will  see  below  that  it  is  most  effective 
as  part  of  a  multiple  time  step  solution  scheme.  For  this  reason  we  will  refer 
to  the  initial  time  as  tk  and  the  hnal  time  as  tk+i  =  tk  +  t.  Second,  the  total 
error  of  the  current  approach  is  separated  into  two  components,  e  =  1  —  6'y, 
where  both  7  and  S  are  numbers  slightly  less  than  1  and  will  be  considered  as 
independent  inputs  to  the  algorithm.  Here  7  refers  to  the  required  sum  of  the 
truncated  probability  distribution  at  tk,  and  S  refers  to  the  relative  accuracy 
requirement  for  the  solution  at  tk+i  compared  to  the  accuracy  at  tk-  Third,  for 
added  convenience  we  will  use  the  notation  E*  =  Vj-  |exp(Aj^r)ej^|  to  de¬ 
note  the  Ji  indexed  FSP  approximation  of  the  distribution  at  tk+i  conditioned 
upon  the  conhguration  at  tk-  Any  matrix  exponential,  exp(A j.r)  provides 
not  only  Ej  but  also  approximations  to  Ej  for  every  j  G  Ji-  Once  we  com¬ 
pute  these  matrix  exponentials,  we  will  store  every  E^  =  T>j.  |exp(Aj-r)e-^, | 
and  its  corresponding  index  set  Jj  =  Ji  that  meets  our  accuracy  requirement 
||Ej||^  >  5-  These  are  then  reused  to  reduce  the  total  number  of  matrix  com¬ 
putations.  With  this  notation,  we  can  now  state  the  following  algorithm; 


The  FSP  Algorithm  for  Non-Sparse  Initial  PDV’s 


Inputs  Propensity  functions  and  stoichiometry  for  all  reactions. 

Error  Parameters,  0  <  7  <  1  and  0  <  5  <  1  . 

Initial  probability  distribution,  P(tfc),  where  1  >  ||P(tfc)||^  >  7 
Length  of  time  interval,  r. 


Step  0  Choose  a  finite  set  of  states,  X4  such  that  ||P4(0)||^  >  7. 
Initialize  a  counter,  i  as  the  first  element  in  Ik- 
Initialize  the  FSP  solution  index  set:  If  =  {i}. 

Initialize  the  FSP  solution  summation  to  zero:  Pf^^^{tf)  =  0. 


Step  1  If  Ei  has  not  already  been  calculated: 

Use  original  FSP  algorithm  to  hnd  Ji  and  exp(Aj.r)  such 


that 


exp(Aj,r)e 


J^ 


For  every  j  &  Ji,  if 


>  5. 


exp(Aji/)ej,  >  5,  then  record 


Ej  =  Vj^  {exp(Aj,t/)e^jJ  and  Jj  =  -Ji 


Step  2  Update  the  FSP  solution  index  set:  If  =  If[jJi- 

Update  the  FSP  solution  summation;  +  PiEi- 


Step  3  If  i  is  the  last  element  in  Jq,  Stop. 

Vj^.  approximates  P(t/)  to  within  e  =  1  —  ■j6. 

Step  4  Increment  i  to  the  next  element  in  Jq  and  return  to  Step  1. 


These  alterations  in  the  FSP  algorithm  enable  one  to  handle  problems  in  which 
the  initial  probability  density  vector  is  not  sparse.  On  its  own,  this  may  be 
convenient  when  we  wish  to  study  systems  that  begin  somewhere  within  a 
range  of  possible  initial  conhgurations.  However,  as  we  will  illustrate  in  the 
following  section,  the  non-sparse  FSP  algorithm  has  its  greatest  use  when  it 
is  integrated  into  a  multiple  time  stepping  FSP  algorithm. 


4  The  Multiple  Time  Step  FSP  Method 


Suppose  that  we  require  that  the  FSP  solution  is  precise  to  a  1-norm  error 
of  e  for  the  entire  time  interval  (0,t/).  This  requires  that  the  system  remains 
with  probability  (1-e)  within  the  hnite  set  Xj  for  all  times  t  G  (0,  tf).  One  can 
envision  many  simple  cases  where  such  a  restriction  can  require  an  exorbitantly 
large  space  Xj.  Suppose  our  system  begins  with  an  initial  condition  at  t  =  0 
far  from  the  support  of  the  distribution  at  the  later  time  as  illustrated  in 
Figure  la.  In  this  case  the  probability  distribution  is  likely  to  evolve  along 
some  path  connecting  the  initial  condition  to  the  hnal  solution.  To  achieve 
acceptable  accuracy  at  all  times,  the  projection  region  must  contain  not  only 
the  initial  condition  and  the  final  solution,  but  also  every  point  likely  to  be 
reached  during  the  intervening  time.  In  such  a  circumstance,  it  can  help  to 
break  the  time  interval  into  pieces  and  require  only  that  the  FSP  criteria  are 
satished  only  during  each  sub-interval.  In  effect,  we  seek  a  changing  projection 
space  that  follows  the  support  of  the  distribution  as  it  evolves.  To  do  this,  we 
utilize  the  linearity  and  time  invariance  properties  of  the  chemical  master 
equation. 


Suppose  we  start  with  a  known  initial  probability  distribution,  P(0),  and  we 
wish  to  approximate  the  solution  to  the  CME  in  k  time  steps  of  equal  length 
r.  Using  the  algorithm  in  the  previous  section,  we  can  specify  a  positive  <5  <  1 
and  require  that  transition  vectors  {Ej}  satisfy  ||Ej||^  >  S  for  all  i.  For  the  first 
time  step,  suppose  that  we  simply  specify  7i  =  5  and  we  use  the  non-sparse 
FSP  algorithm  to  hnd  an  approximation  of  the  distribution  at  =  r  such  that 


0  <  Vj,  <  P(ti)  and  ||pf^^^(U)||^  >  7p5  =  ^2. 

For  the  second  time  step,  we  use  as  the  initial  distribution.  If  we 

use  the  same  6,  we  can  save  some  effort  by  reusing  some  of  the  Ej’s  already 
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computed.  However,  since  our  solution  at  the  end  of  the  previous  step  has  a 
guaranteed  sum  of  only  <5^,  we  must  choose  a  different  72.  A  very  reasonable 
choice  is  simply  to  use  the  guarantee  from  the  previous  step;  72  =  With 
this  choice,  we  can  again  apply  the  non-sparse  FSP  algorithm  to  hnd  an  FSP 
solution  at  the  end  of  the  second  time  step  such  that 


0  <  Vj,  {Pf/^(t2)}  <  P(t2)  and 


> 


Following  this  example,  at  each  step,  if  we  use  7?;  =  <5^,  then  we  will  recover 
a  solution  such  that 


yFSP, 


0  <  »4  {Pf '”(*1)}  <  P(*t) 

If  we  apply  the  fact  that  ||P(ffc)||i  =  1,  we  have 

|PP'’W||,  >(l|Pfe)ll,-l)  +  '5*+'. 

which  after  some  rearranging  yields 

P(4)-I>4{P77r)}||^<l-i‘+'. 


>  5 


fc+l 


Suppose  that  we  wish  to  hnd  a  solution  that  is  within  e  of  the  exact  solution  of 
the  CME  at  tf  =  Kr.  Following  the  ideas  above,  we  would  choose  5  according 
to  the  relation  e  =  1  —  or  <5  =  (1  —  .  This  procedure  is  stated  more 

formally  in  the  following  algorithm. 

The  Multiple  Time  Step  FSP  Algorithm 

Inputs  Propensity  functions  and  stoichiometry  for  all  reactions. 

Initial  probability  distribution,  P(to)- 
Final  time  of  interest,  tf. 

Total  error,  e  >  0. 

Step  0  Choose  the  number  of  time  steps,  iP,  and  calculate  r  =  tf/K. 
Compute  the  required  sum  for  each  E*,  5  =  (1  —  . 

Initialize  time  step  counter:  k=0. 

Choose  initial  time  index  set,  Jq,  such  that  ||P/o(^o)|li  > 

Initialize  the  FSP  approximate  solution  at  to  ,Pfo^^(to)=P/o(M- 

Step  1  Run  the  Non-Sparse  FSP  algorithm  with  the  initial  condition 
and  error  parameters  5  and  7^  =  and  get  Pf^_^^(tk+i). 

Step  2  If  k  +  1  =  K,  then  Stop. 

77 approximates  Pi^,{tf)  to  within  e. 
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(C)  (d) 

Fig.  1.  Schematic  of  the  Multiple  Time  Step  FSP  method,  (a)  We  are  given  a 
Markov  process  that  begins  at  a  known  initial  point  in  the  configuration  space. 
As  the  probability  distribution  evolves,  it  follows  a  long  path  in  the  configuration 
space  such  that  at  time  te  the  distribution  is  supported  in  a  region  far  from  the 
initial  condition,  (b)  In  order  to  find  a  sufficiently  accurate  FSP  solution  for  all 
times  in  the  interval  [0, 6r],  the  FSP  must  include  not  only  the  initial  condition  and 
the  final  distribution,  but  also  all  points  along  the  path,  (c)  To  save  computational 
effort,  one  can  discretize  the  time  interval  into  smaller  intervals  and  find  overlapping 
projections  that  need  only  satisfy  the  accuracy  requirements  during  those  shorter 
periods  of  time.  Here  the  final  distribution  of  each  time  step  (shown  in  grey)  becomes 
the  initial  distribution  for  the  next  time  step  (shown  in  black),  (d)  The  end  result 
is  a  discrete  map  taking  the  distribution  from  one  instant  in  time  to  the  next. 
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Step  3  Increment  k  and  return  to  Step  1. 


To  see  how  one  may  benefit  from  this  modihcation  to  the  FSP  algorithm,  we 
refer  to  Figure  1.  Suppose  that  we  are  interested  in  hnding  the  distribution 
at  time  t  =  6r  of  a  Markov  process  that  begins  in  the  known  initial  config¬ 
uration  represented  by  the  black  dot.  Even  though  the  distribution  at  each 
of  the  times  {0,  r,  2r, . . . ,  6r}  are  supported  on  only  a  small  portion  of  the 
conhguration  space,  the  one  shot  FSP  solution  must  include  the  whole  region 
of  the  conhguration  space  that  is  swept  by  the  distribution  between  0  and  6r 
(see  Fig  lb).  Therefore,  the  one  step  FSP  algorithm  requires  a  large  matrix 
exponential  computation.  By  splitting  the  full  interval  into  six  subintervals  as 
shown  in  Fig.  Ic,  we  will  require  more  exponential  computations,  but  since 
each  of  these  computations  will  be  much  smaller,  the  total  computational  ef¬ 
fort  may  be  much  less.  The  following  subsection  illustrates  the  use  of  this 
algorithm  through  a  simplihed  model  of  the  heat  shock  response  in  E-coli. 


4-1  Toy  Heat  Shock  Example 


When  a  cell’s  environment  changes,  that  cell  must  either  adapt  or  perish. 
Continually  faced  with  this  choice,  life  has  evolved  to  contain  many  complex 
gene  regulatory  networks  that  allow  for  quick  and  precise  adaptation.  The 
heat  shock  response  in  E.  coli  is  excellent  example  of  one  such  mechanism 
[20].  A  simplihed  version  of  this  system  consists  of  three  biochemical  species 
that  interact  according  to  a  set  of  three  reactions. 


-Si  ^  S2  ^  S3,  (18) 

where  si,  S2  and  S3  correspond  to  the  (T32-DnaK  complex,  the  (T32  heat  shock 
regulator  and  the  (T32-RNAP  complex,  respectively.  This  model  of  the  heat 
shock  subsystem  has  been  analyzed  before  using  various  computational  meth¬ 
ods  including  Monte  Carlo  implementations  [7,21]  as  well  as  the  FSP  with  the 
multiple  time-scale  model  reduction  [17].  Here  we  use  this  model  in  order  to 
illustrate  our  current  computational  algorithm. 

We  will  take  the  reaction  rates  and  initial  conhguration  for  the  system  to  be: 

Cl  =  10,  C2  =  4  X  10^  C3  =  2,  si(0)  =  2000,  S2(0)  =  S3(0)  =  0.  (19) 

For  these  parameters  there  are  2,001,000  reachable  conhgurations,  and  the  full 
chemical  master  equation  is  far  too  large  to  be  solved  exactly.  Applying  the 
Finite  State  Projection  method  allows  us  to  signihcantly  reduce  the  order  of 
the  problem  and  achieve  a  manageable  solution  at  least  for  small  time  intervals 
{t  <  300s).  In  order  to  retain  accuracy  in  the  solution  to  an  error  of  e  =  10“^ 
at  a  time  tf  =  100s,  we  need  include  the  1430  conhgurations  where  S3  <  129 
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Fig.  2.  Probability  distribution  of  the  population  of  the  1T32-RNAP  complex  at  100, 
200  and  300  seconds  as  computed  with  the  basic  FSP  algorithm  (solid  line)  and 
with  the  Multiple  Time  Step  FSP  algorithm  (dotted  line). 

and  S2  <  10.  For  larger  time  intervals,  we  need  to  inclnde  more  points  in 
the  FSP  solntion:  for  t/  =  200s,  we  need  inclnde  2585  confignrations  where 
S3  <  234  and  S2  <  10;  and  for  tf  =  300s,  we  need  inclnde  3641  conhgnrations 
where  S3  <  330  and  S2  <  10.  The  solid  lines  in  Fignre  2  present  the  probability 
distribntions  for  the  nnmber  of  S3  molecnles  at  100,  200  and  300s,  and  the  top 
section  of  Table  1  snmmarizes  the  compntational  reqnirements  necessary  to 
achieve  these  resnlts  with  the  basic  FSP  algorithm.  We  have  also  applied  the 
mnltiple  time  step  FSP  algorithm  to  the  toy  heat  shock  model  in  order  to 
compnte  the  probability  distribntion  at  these  same  times.  As  a  first  attempt, 
we  have  nsed  a  step  size  of  five  seconds.  Fignre  2  shows  that  there  is  no 
discernible  difference  between  the  resnlts  of  the  basic  FSP  algorithm  and  those 
of  the  mnltiple  time  step  FSP  algorithm.  For  a  comparison  of  compntational 
efforts,  the  bottom  half  of  Table  1  provides  the  maximnm  matrix  size,  nnmber 
of  matrices,  and  compntational  time  reqnired  for  the  mnltiple  time  step  FSP 
algorithm.  As  the  total  time  increases  from  100  to  300  seconds,  so  does  the 
compntational  benefit  of  the  discretized  algorithm;  for  a  final  time  of  100,  the 
cnrrent  algorithm  rednces  compntational  time  by  a  factor  of  abont  2.5,  for  a 
hnal  time  of  300,  the  redaction  is  a  factor  of  abont  11.8. 

Of  conrse,  while  the  accnracy  of  the  mnltiple  time  step  FSP  is  gnaranteed,  the 
efficiency  of  the  algorithm  will  nndonbtedly  depend  npon  onr  chosen  step  size. 
Fignre  3  illnstrates  some  of  the  snbtleties  of  this  tradeoff  by  plotting  the  size  of 
the  largest  exponentiated  matrix,  the  nnmber  of  matrix  exponentials,  and  the 
compntational  time  all  as  fnnctions  of  the  nnmber  of  time  steps  (bottom  axis) 
and  the  step  length  (top  axis).  As  we  nse  more  time  steps,  the  probability 
distribntion  has  less  time  to  disperse  between  one  step  and  the  next,  and 
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the  required  matrix  exponentials  are  smaller  as  shown  in  Fig.  3a.  However, 
because  the  matrix  dimension  is  a  discrete  integer  quantity,  this  decrease  is 
stepwise  rather  than  smooth,  and  a  large  range  of  step  lengths  may  require 
the  same  matrix  size.  If  a  step  length  is  at  the  low  end  of  that  range,  the 
matrix  exponentials  required  to  get  each  E*  are  often  slightly  more  precise 
than  is  absolutely  necessary,  and  are  therefore  more  likely  to  provide  other 
Ej’s  as  well-fewer  exponential  computations  are  necessary.  Conversely,  if  a 
step  length  is  at  the  high  end  of  the  range  for  a  given  matrix  size,  fewer  Ej’s 
will  come  from  each  exponential  computation-more  exponential  computations 
are  necessary.  This  trend  is  clear  when  one  compares  Fig  3a  to  3b. 

In  order  to  show  how  these  concerns  affect  the  computation,  we  have  broken 
the  total  computational  cost  in  Fig  3c;  into  two  components.  The  hrst  cost  is 
that  of  computing  the  matrix  exponentials.  The  second  cost  is  that  of  updating 
the  solution  from  one  step  to  the  next.  For  tj  =  300s,  this  tradeoff  is  optimized 
for  130  time  steps  corresponding  to  a  interval  length  of  r  ~  2.3s.  To  obtain 
the  solution  with  this  time  step,  the  algorithm  needed  to  compute  50  matrix 
exponentials  of  size  198  x  198  or  smaller,  and  the  computation  takes  about 
28s. 

Extrapolating  from  Table  l(top)  suggests  that  a  regular  FSP  solution  at 
t  =  1000s  would  require  inclusion  of  more  than  11000  conhgurations.  Un¬ 
fortunately,  the  memory  requirement  to  exponentiate  a  11000  x  11000  matrix 
exceeds  the  specihcations  of  our  machine,  and  the  unmodihed  FSP  algorithm 
cannot  be  used.  Alternatively,  by  discretizing  the  full  time  interval,  we  can 
signihcantly  reduce  the  computational  complexity  and  bring  the  model  back 
into  the  realm  of  a  solvable  problem.  Once  again,  there  is  a  dehnite  tradeoff 
between  too  many  and  too  few  time  steps,  and  Figure  4  plots  the  the  size  of 
the  largest  exponentiated  matrix,  the  number  of  matrix  exponentials,  and  the 
computational  time  as  a  function  of  the  number  of  time  steps.  For  tf  =  1000s, 
the  computational  tradeoff  is  optimized  for  290  time  steps  corresponding  to  a 
interval  length  r  ~  3.4.  To  obtain  the  solution  with  this  time  step,  the  algo¬ 
rithm  needed  to  compute  107  matrix  exponentials  of  size  252  x  252  or  smaller, 
and  the  computation  takes  about  166s. 


5  Conclusions 


Although  the  original  hnite  state  projection  method  can  signihcantly  reduce 
the  order  of  the  chemical  master  equation  for  many  problems,  this  initial  re¬ 
duction  is  not  sufficient  for  all  systems.  Fortnnately,  the  FSP  is  amenable  to 
numerous  modihcations,  which  can  considerably  improve  upon  the  method’s 
range  and  potency.  In  this  paper  we  have  concentrated  on  one  computational 
difficnlty  that  arises  when  system  trajectories  slowly  drift  over  large  regions  of 
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Fig.  3.  Trade  off  between  more  and  fewer  time  steps  in  the  Multiple  Time  Step 
FSP  (MTS-FSP)  algorithm  solution  for  the  toy  heat  shock  model  at  a  final  time  of 
tf  =  300s.  The  following  are  plotted  as  function  of  the  number  of  steps:  (top)  the 
size  of  the  largest  required  matrix  exponential  computation,  (middle)  the  number 
of  matrix  exponential  computations  performed,  (bottom)  the  computational  time 
required  for  the  MTS-FSP  algorithm  split  into  two  components:  the  exponential 
computation  costs,  and  the  overhead  costs.  All  computations  have  been  performed 
in  Matlab  7.2  on  a  Dual  2  Ghz  PowerPC  G5. 
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Fig.  4.  Trade  off  between  more  and  fewer  time  steps  in  the  Multiple  Time  Step 
FSP  (MTS-FSP)  algorithm  solution  for  the  toy  heat  shock  model  at  a  final  time  of 
tf  =  1000s.  The  following  are  plotted  as  function  of  the  number  of  steps:  (top)  the 
size  of  the  largest  required  matrix  exponential  computation,  (middle)  the  number 
of  matrix  exponential  computations  performed,  (bottom)  the  computational  time 
required  for  the  MTS-FSP  algorithm  split  into  two  components:  the  exponential 
computation  costs,  and  the  overhead  costs.  All  computations  have  been  performed 
in  Matlab  7.2  on  a  Dual  2  Ghz  PowerPC  G5. 
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Table  1 

Comparison  of  the  computational  requirements  of  the  basic  FSP  algorithm  and  the 
Multiple  Time  Step  FSP  (MTS-FSP)  algorithm  for  the  toy  heat  shock  example 
for  three  different  final  times:  tf  =  100,  200,  and  300  seconds.  For  the  MTS-FSP 
algorithm  results  in  this  table,  we  have  used  a  non-optimal  step  size  of  5  seconds, 
different  step  sizes  will  provide  the  same  accuracy  at  a  lower  computational  cost. 


Original  FSP  Algorithm 

Final  Time  (s) 

Matrix  Size 

Comp.  Time®'  (s) 

Error  .  ^ 

100 

1430  X  1430 

27 

9.61  X  10“^ 

200 

2585  X  2585 

165 

9.04  X  10“^ 

300 

3641  X  3641 

437 

9.67  X  10“^ 

1000 

11000  X  11000 

Multiple  Time  Step  FSP  Algorithm 

Final  Time  (s) 

^  Exponentials 

Matrix  Size 

Comp.  Time  (s) 

Error  . 

100 

14 

275  X  275 

11 

2.84  X  10“^ 

200 

27 

275  X  275 

23 

4.25  X  10“^ 

300 

39 

275  X  275 

37 

5.08  X  10“^ 

1000 

104 

252  X  252 

180 

1.8  X  10“*^ 

^  Computations  have  been  performed  in  Mat- 
lab  7.2  on  a  Dual  2  GHz  PowerePC  G5 


the  configuration  space  during  long  time  intervals.  In  order  to  use  the  original 
FSP  method  for  these  cases,  one  must  include  vast  portions  of  the  conhgura- 
tion  space  in  the  projected  solution.  As  the  size  of  the  included  conhguration 
space  increases,  so  do  the  computational  requirements  of  the  FSP.  However, 
in  some  cases  this  difficulty  can  be  ameliorated  simply  by  solving  the  FSP  for 
a  series  of  smaller  time  intervals.  Here  we  have  presented  the  Multiple  Time 
Step  FSP  (MTS-FSP)  algorithm,  which  is  essentially  an  incremental  approach 
to  solving  the  original  FSP. 

The  MTS-FSP  algorithm  is  built  upon  three  important  aspects  that  the  FSP 
inherits  from  the  chemical  master  equation;  linearity,  time-invariance,  and  pos¬ 
itivity.  The  linearity  of  the  FSP  allows  us  to  apply  the  principle  of  superposi¬ 
tion  with  regards  to  initial  conditions-if  we  know  the  probability  distribution 
at  time  0  and  we  know  the  conditional  probabilities  at  time  r  conditioned 
on  each  conhguration  at  time  0,  then  we  can  easily  compute  the  probability 
distribution  at  time  r.  The  time  invariance  of  the  FSP  assures  us  that  if  we 
know  the  probabilities  at  time  r  conditioned  on  time  0,  then  we  also  know  the 
probabilities  at  time  t  +  r  conditioned  on  any  time  t.  The  positivity  of  the  FSP 
guarantees  us  that  we  never  over-predict  the  solution  to  the  CME.  Whether 
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we  neglect  some  portion  of  the  initial  probability  distribution  or  lose  some  of 
that  distribution  to  conhgurations  excluded  from  our  various  projections,  the 
resulting  error  is  known  at  the  end  of  each  time  step.  By  directly  controlling 
the  error  at  each  time-step,  we  can  control  the  hnal  error. 

We  have  demonstrated  the  MTS-FSP  algorithm  on  the  toy  heat  shock  prob¬ 
lem.  For  time  intervals  of  one,  two  and  three  hundred  seconds  the  FSP  and  the 
current  algorithm  produce  nearly  identical  results,  but  with  the  new  method, 
we  can  compute  those  results  much  faster.  In  addition,  the  new  algorithm 
extends  the  range  of  problems  to  which  the  FSP  approach  may  be  applied. 
To  solve  the  toy  heat  shock  problem  over  a  time  interval  of  one  thousand 
seconds,  original  FSP  algorithm  requires  a  conhguration  space  of  over  11000 
points,  which  is  too  large  to  manage.  However,  we  can  now  easily  solves  the 
problem  using  matrices  less  than  one  thirtieth  of  the  size. 

This  time  stepping  approach  is  just  one  of  many  mutually  benehcial  improve¬ 
ments  that  are  quickly  expanding  the  ability  of  the  FSP  to  directly  solve  the 
chemical  master  equation.  This  current  approach  retains  the  full  accuracy  and 
properties  of  the  original  FSP  and  can  easily  be  combined  with  other  model 
reduction  techniques  such  as  those  based  linear  systems  and  modern  control 
theory.  While  we  may  never  be  able  to  directly  solve  every  master  equation, 
it  remains  to  be  seen  just  how  far  these  FSP  based  approaches  can  push  back 
the  boundary  between  solvable  and  unsolvable. 
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